Última actualización: 25 mayo 2020





Generalidades

En esta entrega se describe la estrategia de calibración usada con 5 estaciones de caudal. Los datos observados de caudal han sido previamente tratados (ver: Tratamiento de caudales).

Periodo de calibración y validación

La siguiente tabla muestra los periodos de validación.

t23 2004-01-01 2005-12-31
t30 2005-01-01 2006-12-31
t36 2008-01-01 2009-12-31
t38 2007-01-01 2008-12-31
t41 2006-01-01 2007-12-31

Las siguientes figuras muestran en gris los periodos de validación, el resto de la serie se utiliza para calibración.

fig = plotly_empty()

fig = fig %>% add_segments(x = ~t23[1], y=max(qobs$q23,na.rm = T), yend = max(qobs$q23,na.rm = T), 
                           xend = t23[2], name="Validación Paso Pache",
                           line = list(color = 'black', width = 2))

fig = fig %>% add_segments(x = ~t30[1], y=max(qobs$q30,na.rm = T), yend = max(qobs$q30,na.rm = T), 
                           xend = t30[2], name="Validación Paso Roldán",
                           line = list(color = 'blue', width = 2))

fig = fig %>% add_segments(x = ~t36[1], y=max(qobs$q36,na.rm = T), yend = max(qobs$q36,na.rm = T), 
                           xend = t36[2], name="Validación Fray Marcos",
                           line = list(color = 'red', width = 2))

fig = fig %>% add_segments(x = ~t38[1], y=max(qobs$q38,na.rm = T), yend = max(qobs$q38,na.rm = T),
                           xend = t38[2], name="Validación San Ramon",
                           line = list(color = 'orange', width = 2))

fig = fig %>% add_segments(x = ~t41[1], y=max(qobs$q41,na.rm = T), yend = max(qobs$q41,na.rm = T),
                           xend = t41[2], name="Validación Paso Pache",
                           line = list(color = 'green', width = 2))

fig <- fig %>%add_lines(x = ~qobs$Rdate, y = ~qobs$q23, type="scatter",mode = "lines",
               name="Paso de los Troncos", line = list(color = 'black', width = 2))

  
fig <- fig %>%add_lines(x = ~qobs$Rdate, y = ~qobs$q30, type="scatter",mode = "lines",
               name="Paso Roldán", line = list(color = 'blue', width = 2))

  
fig <- fig %>%add_lines(x = ~qobs$Rdate, y = ~qobs$q36, type="scatter",mode = "lines",
               name="Fray Marcos", line = list(color = 'red', width = 2))

  
fig <- fig %>%add_lines(x = ~qobs$Rdate, y = ~qobs$q38, type="scatter",mode = "lines",
               name="Paso San Ramón", line = list(color = 'orange', width = 2))

  
fig <- fig %>%add_lines(x = ~qobs$Rdate, y = ~qobs$q41, type="scatter",mode = "lines",
               name="Paso Pache", line = list(color = 'green', width = 2))



fig <- fig %>% layout(#legend = list(x = 1, y = 0.9),
                      yaxis = list(title="Runoff (mm)",
                                   showticklabels = TRUE, ticks="outside",
                                   showline=TRUE, showgrid=TRUE),
                      xaxis = list(title="",showticklabels = TRUE, 
                                   ticks="outside",showline=TRUE, 
                                   showgrid=TRUE, range=as.Date(c("2004-01-01", "2012-12-31"))),
                      height=300,
                      width=700)

  
fig

Anexo

rm(list = ls())

# Librerias
library(SWATplusR)
library(stringr)
library(hydroGOF)
library(hydromad)
library(tibble)
library(hydroPSO)

# Caudales observados
obs.data <- readRDS("./data/qobs_5est_2004_2012_para_calibracion.RDS")

# ajuste de caudal
obs.data$q23 = 1.630*obs.data$q23
obs.data$q30 = 0.789*obs.data$q30
obs.data$q36 = 1.184*obs.data$q36
obs.data$q38 = 0.724*obs.data$q38
obs.data$q41 = 1.120*obs.data$q41


# SWAT project mensual
swat_project = "./SWAT_SL_20200401/"

# Funcion de optimización
swat_optim <- function(mypar, obs.data, project_dir, start_sim, end_sim,runrun_path, head_log) {
  
  
  names(mypar) = head_log[25:length(head_log)]
  path_run = paste0(runrun_path,"run_",Sys.getpid())
  
  par_set = tibble("CN2::CN2.mgt|change = pctchg" = mypar["CN2"]*100,
                   
                   "ESCO::ESCO.hru|change = absval" = mypar["ESCO"],
                   "SOL_AWC::SOL_AWC.sol|change = pctchg " = mypar["SOL_AWC"]*100,
                   "ALPHA_BF::ALPHA_BF.gw|change = absval" = mypar["ALPHA_BF"],
                   "GWQMN::GWQMN.gw|change = absval" = mypar["GWQMN"],
                   "GW_DELAY::GW_DELAY.gw|change = absval |" = mypar["GW_DELAY"],
                   "REVAPMN::REVAPMN.gw|change = absval" =  mypar["REVAPMN"],
                   "GW_REVAP::GW_REVAP.gw|change = absval" = mypar["GW_REVAP"],
                   "OV_N::OV_N.hru|change = absval" = mypar["OV_N"],
                   "SLSUBBSN::SLSUBBSN.hru|change = pctchg" = mypar["SLSUBBSN"]*100,
                   
                   
                   "HRU_SLP::HRU_SLP.hru|change= pctchg" = mypar["HRU_SLP"]*100,
                   "USLE_K::USLE_K.sol|change= pctchg" = mypar["USLE_K"]*100,
                   "USLE_P::USLE_P.mgt|change= pctchg"=mypar["USLE_P"]*100,
                   "P_UPDIS::P_UPDIS.bsn|change= absval"= mypar["P_UPDIS"], 
                   "FRT_SURFACE::FRT_SURFACE.mgt|change=absval" = mypar["FRT_SURFACE"], 
                   "FRT_KG::FRT_KG.mgt|change=pctchg" = mypar["FRT_KG"]*100,
                   "NPERCO::NPERCO.bsn|change= absval"=mypar["NPERCO"],
                   "ERORGN::ERORGN.hru|change= absval"= mypar["ERORGN"], 
                   
                   "PPERCO::PPERCO.bsn|change= absval"= mypar["PPERCO"],                    
                   "CDN::CDN.bsn|change= absval"= mypar["CDN"],                    
                   "CMN::CMN.bsn|change= absval"= mypar["CMN"],                    
                   "PSP::PSP.bsn|change= absval"= mypar["PSP"]
  ) 
  
  
  # run swat
  qsim = run_swat2012(project_path = project_dir,
                      output = list(FLOW = define_output(file = "rch",
                                                         variable = "FLOW_OUT",
                                                         unit = c(23,30,36,38,41)),
                                    PT = define_output(file = "rch",
                                                       variable = "TOT_Pkg",
                                                       unit = c(23,41)),
                                    NO3 = define_output(file = "rch",
                                                        variable = "NO3_OUT",
                                                        unit = c(23,41))),
                      parameter = par_set,
                      start_date = as.Date("2000-01-01"),
                      end_date = as.Date("2012-12-31"),
                      years_skip = 4,
                      run_path=path_run,
                      quiet = TRUE,
                      output_interval="d",
                      add_parameter = FALSE)
  
  # objective function
  nse23 = NSE(as.vector(qsim$FLOW_23), as.vector(obs.data$q23))
  nse30 = NSE(as.vector(qsim$FLOW_30), as.vector(obs.data$q30))
  nse36 = NSE(as.vector(qsim$FLOW_36), as.vector(obs.data$q36))
  nse38 = NSE(as.vector(qsim$FLOW_38), as.vector(obs.data$q38))
  nse41 = NSE(as.vector(qsim$FLOW_41), as.vector(obs.data$q41))
  
  kge23 = KGE(as.vector(qsim$FLOW_23), as.vector(obs.data$q23))
  kge30 = KGE(as.vector(qsim$FLOW_30), as.vector(obs.data$q30))
  kge36 = KGE(as.vector(qsim$FLOW_36), as.vector(obs.data$q36))
  kge38 = KGE(as.vector(qsim$FLOW_38), as.vector(obs.data$q38))
  kge41 = KGE(as.vector(qsim$FLOW_41), as.vector(obs.data$q41))
  
  pbias23 = pbias(as.vector(qsim$FLOW_23), as.vector(obs.data$q23))
  pbias30 = pbias(as.vector(qsim$FLOW_30), as.vector(obs.data$q30))
  pbias36 = pbias(as.vector(qsim$FLOW_36), as.vector(obs.data$q36))
  pbias38 = pbias(as.vector(qsim$FLOW_38), as.vector(obs.data$q38))
  pbias41 = pbias(as.vector(qsim$FLOW_41), as.vector(obs.data$q41))
  
  
  # objective quality
  qlt_41 = data.frame(flow_lps_km2=1000*qsim$FLOW_41/4896,
                      PT_mgs_km2 = 1000000*qsim$PT_41/(24*60*60*4896),
                      NO3_mgs_km2 = 1000000*qsim$NO3_41/(24*60*60*4896),
                      qobs_lps_km2 = 1000*obs.data$q41/4896,
                      PTobs_mgs_km2 = 0.16*(1000*obs.data$q41/4896)^1.02,
                      NO3obs_mgs_km2 = 0.1*(1000*obs.data$q41/4896)^1.41)
  
  qlt_41[qlt_41<=0] = NA
  qlt_41 = na.omit(qlt_41)
  
  qlt_23 = data.frame(flow_lps_km2=1000*qsim$FLOW_23/687,
                      PT_mgs_km2 = 1000000*qsim$PT_23/(24*60*60*687),
                      NO3_mgs_km2 = 1000000*qsim$NO3_23/(24*60*60*687),
                      qobs_lps_km2 = 1000*obs.data$q23/687,
                      PTobs_mgs_km2 = 0.06*(1000*obs.data$q23/687)^0.93,
                      NO3obs_mgs_km2 = 0.04*(1000*obs.data$q23/687)^1.41)
  
  qlt_23[qlt_23<=0] = NA
  qlt_23 = na.omit(qlt_23)
  
  fit = lm(log(PT_mgs_km2)~log(flow_lps_km2), data=qlt_41)
  a = exp(fit$coefficients[1])
  b = fit$coefficients[2]
  qlt_41$PT_fit = a*(qlt_41$flow_lps_km2)^b
  
  fit = lm(log(NO3_mgs_km2)~log(flow_lps_km2), data=qlt_41)
  a = exp(fit$coefficients[1])
  b = fit$coefficients[2]
  qlt_41$NO3_fit = a*(qlt_41$flow_lps_km2)^b
  
  fit = lm(log(PT_mgs_km2)~log(flow_lps_km2), data=qlt_23)
  a = exp(fit$coefficients[1])
  b = fit$coefficients[2]
  qlt_23$PT_fit = a*(qlt_23$flow_lps_km2)^b
  
  fit = lm(log(NO3_mgs_km2)~log(flow_lps_km2), data=qlt_23)
  a = exp(fit$coefficients[1])
  b = fit$coefficients[2]
  qlt_23$NO3_fit = a*(qlt_23$flow_lps_km2)^b
  
  kge_PT_41 = KGE(qlt_41$PT_fit, qlt_41$PTobs_mgs_km2)
  kge_NO3_41 = KGE(qlt_41$NO3_fit, qlt_41$NO3obs_mgs_km2)
  kge_PT_23 = KGE(qlt_23$PT_fit, qlt_23$PTobs_mgs_km2)
  kge_NO3_23 = KGE(qlt_23$NO3_fit, qlt_23$NO3obs_mgs_km2)
  
  pbias_PT_23 = pbias(qlt_23$PT_fit, qlt_23$PTobs_mgs_km2)
  pbias_PT_41 = pbias(qlt_41$PT_fit, qlt_41$PTobs_mgs_km2)
  
  pbias_NO3_23 = pbias(qlt_23$NO3_fit, qlt_23$NO3obs_mgs_km2)
  pbias_NO3_41 = pbias(qlt_41$NO3_fit, qlt_41$NO3obs_mgs_km2)
  
  
  FO = 0.5*(kge23 + kge30 + kge36 + kge38 + kge41)/5 + 0.25*(kge_PT_41 + kge_PT_23)/2 + 0.25*(kge_NO3_41 + kge_NO3_23)/2 
  
  # logfile
  vout = data.frame(FO,kge23,kge30,kge36,kge38,kge41,
                    nse23,nse30,nse36,nse38,nse41,
                    pbias23,pbias30,pbias36,pbias38,pbias41,
                    kge_PT_23, kge_PT_41, kge_NO3_41, kge_NO3_23,
                    pbias_PT_23, pbias_PT_41, pbias_NO3_23, pbias_NO3_41,
                    t(mypar))
  colnames(vout) = head_log
  log.file = paste0(runrun_path,"log_",Sys.getpid(),".txt")
  
  if(!any(list.files(runrun_path)==paste0("log_",Sys.getpid(),".txt"))){
    write.table(t(head_log), log.file, col.names = FALSE, row.names = FALSE, quote = FALSE)
  }
  
  ff = read.table(log.file, col.names = head_log,colClasses = "character")
  ff = rbind(ff,vout)
  write.table(ff, log.file, col.names = FALSE, row.names = FALSE, quote = FALSE)
  
  #output
  return(FO)
}

CN2 = c(-0.2,0.2)

ESCO = c(0.1,0.45)
SOL_AWC = c(-0.2,.2)
ALPHA_BF = c(0.4,0.8)  
GWQMN = c(800,4000)
GW_DELAY = c(20,45)
REVAPMN = c(300,700)
GW_REVAP = c(0.02, 0.01)
OV_N = c(0.4, 0.6)
SLSUBBSN = c(-0.1, 0.3)

HRU_SLP = c(-0.1,0.1)
USLE_K = c(-0.2,0.2)
USLE_P = c(-0.2,0.2)
P_UPDIS = c(20, 90)
FRT_SURFACE = c(0.1,0.6)
FRT_KG = c(-0.3,0.3)
NPERCO = c(0.5,1)
ERORGN = c(2,4)
PPERCO = c(10,17.5) #bsn absval
CDN = c(0,3) #bsn absval
CMN = c(0,0.03) #bsn absval
PSP = c(0.01,0.8)


names_par = c("CN2","ESCO","SOL_AWC",
              "ALPHA_BF","GWQMN","GW_DELAY","REVAPMN","GW_REVAP","OV_N","SLSUBBSN","HRU_SLP",
              "USLE_K","USLE_P","P_UPDIS","FRT_SURFACE","FRT_KG","NPERCO","ERORGN","PPERCO","CDN","CMN","PSP")

#           CN2   ESCO  SOL_AWC ALPHA_BF  GWQMN GW_DELAY  REVAPMN GW_REVAP  OV_N  SLSUBBSN
par_lwr = c(CN2[1],
            ESCO[1], SOL_AWC[1], ALPHA_BF[1], GWQMN[1], GW_DELAY[1], REVAPMN[1], GW_REVAP[1], OV_N[1],SLSUBBSN[1],
            HRU_SLP[1], USLE_K[1], USLE_P[1], P_UPDIS[1], FRT_SURFACE[1], FRT_KG[1], NPERCO[1], ERORGN[1], PPERCO[1], CDN[1], CMN[1],PSP[1])

par_upr = c(CN2[2],
            ESCO[2], SOL_AWC[2], ALPHA_BF[2], GWQMN[2], GW_DELAY[2], REVAPMN[2], GW_REVAP[2], OV_N[2],SLSUBBSN[2],
            HRU_SLP[2], USLE_K[2], USLE_P[2], P_UPDIS[2], FRT_SURFACE[2], FRT_KG[2], NPERCO[2], ERORGN[2], PPERCO[2], CDN[2], CMN[2],PSP[2])

par_init = (par_upr + par_lwr)/2


# Optimización
head_log = c("FO","KGE23","KGE30","KGE36", "KGE38","KGE41",
             "NSE23","NSE30","NSE36","NSE38", "NSE41",
             "pbias23","pbias30","pbias36","pbias38", "pbias41",
             "kge_PT_23", "kge_PT_41", "kge_NO3_41", "kge_NO3_23",
             "pbias_PT_23", "pbias_PT_41", "pbias_NO3_23", "pbias_NO3_41",
             names_par)


fr = list.files("./run_optim_mensual/", pattern = "*.txt", full.names=TRUE)
file.remove(fr)
opt_pso = hydroPSO(par=par_init, fn = swat_optim,lower=par_lwr, upper=par_upr, 
                   control = list(MinMax = "max",
                                  maxit=500,
                                  reltol=0.00001,
                                  parallel="parallelWin", 
                                  par.nnodes=12,
                                  par.pkgs = list("hydromad", 
                                                  "SWATplusR",
                                                  "stringr",
                                                  "hydroGOF",
                                                  "tibble"),
                                  write2disk=FALSE,
                                  REPORT=1,
                                  npart=48,
                                  normalise=TRUE),
                   obs.data = obs.data, 
                   project_dir = swat_project, 
                   runrun_path = "./run_optim_mensual/",
                   head_log = head_log
)
 




A work by Rafael Navas